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Abstract 

Complex systems display variability over a broad range of spatial and temporal 
scales. Some scales are unresolved due to computational limitations. The impact 
of these unresolved scales on the resolved scales needs to be parameterized or taken 
into account. One stochastic parameterization scheme is devised to take the effects 
of unresolved scales into account, in the context of solving a nonlinear partial dif- 
ferential equation with memory (a time- integral term), via large eddy simulations. 
The obtained large eddy simulation model is a stochastic partial differential equation. 
Numerical experiments are performed to compare the solutions of the original system 
and of the stochastic large eddy simulation model. 
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1 Introduction 



Stochastic parameterizations have been investigated intensively for quantifying uncertain- 
ties in mathematical models of physical, geophysical, environmental, and biological sys- 
tems. We may roughly classify such uncertainties into two kinds. The first kind of un- 
certainties may be called model uncertainty. They involve with physical processes that 
are less known, not yet well understood, not well-observed or measured, difficult to be 
described in the mathematical models, or otherwise ignored in the usual deterministic 
modeling. Parameterizations of model uncertainties have been considered in, for example, 
[HI EU [361 IS [Ml E3 [26] and references therein. 

The second kind of uncertainties may be called simulation uncertainty. This arises 
in numerical simulations of multiscale systems that display a wide range of spatial and 
temporal scales, with no clear scale separation. Due to the limitations of computer power, 
at present and for the conceivable future, not all scales of variability can be explicitly 
simulated or resolved. Although these unresolved scales may be very small or very fast, 
their long time impact on the resolved simulation may be delicate (i.e., may be negligible 
or may have significant effects, or in other words, uncertain). Thus, to take the effects of 
unresolved scales on the resolved scales into account, representations or parameterizations 
of these effects are required [30l [3] . 

The present paper deals with simulation uncertainty, i.e., stochastically parameterizing 
the effects of the unresolved scales on the resolved scales. We consider this issue in the 
context of large eddy simulations (LES) of a nonlinear partial differential equation with 
memory. Relevant existing works include (IC1H1H51H21II51I5HI5|I51I57|I55]. 

In large eddy simulations of fluid or geophysical fluid flows [31 [30] , the unresolved 
scales appear as the so-called subgrid scales (SGS). The SGS term appears to be highly 
fluctuating ("random"); see the Figure 1 in [22]. Partially motivated by this, stochastic 
parameterizations of subgrid scales have been investigated in fluid, geophysical and climate 
simulations, based on physical or intuitive or empirical arguments. Another, perhaps more 
important, motivation for applying stochastic parameterizations of subgrid scales is to 
induce the desired backward energy flux ("stochastic backscatter" ) in fluid simulations 

PH OS [32]. 

We present one stochastic parameterization scheme of the subgrid scale term in the 
large eddy simulation of a nonlinear partial differential equation with an extra memory 
term, which is in fact a nonlinear integro-partial differential equation. The approximation 
scheme is based on stochastic calculus involving with a fractional Brownian Motion, and 
the "parameter' to be calculated is a spatial function, which is derived using Ito stochastic 
calculus. 

This paper is organized as follows. After introducing large eddy simulations in §2, 
we discuss fractional Brownian motions and colored noise in §3, and devise a stochastic 
parameterization scheme of the subgrid scales in details in §3 and §4, respectively. Fi- 
nally, in §5, we demonstrate this stochastic parameterization scheme by a few numerical 
experiments on solving a nonlinear partial differential equation with memory. 
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2 Stochastic large eddy simulations 



As an example, we consider the following nonlinear partial differential equation with a 
memory term (time-integral term) 

/"* 1 

ut = u xx + u-u 3 + J Y+T£Z g|/3 U ( X i S ) ds ' (!) 

under appropriate initial condition u(x,0) = uq(x) and boundary conditions u( — l,t) = 
a, u(l,t) = b with a, b constants, on a bounded domain D : — 1 < x < 1. Here (3 is 
a positive constant. This model arises in mathematical modeling in ecology [39], heat 
conduction in certain materials [T7] and materials science [T7j. The time- integral 
term here represents a memory effect depending on the past history of the system state, 
and this memory effect decays polynomially fast in time. 

The large eddy solution u is the true solution u looked through a filter: i.e., through 
convolution with a spatial filter G$(x), with spatial scale (or filter size or cut-off size) 
5 > 0: 

u(x,t) := u* G s = / u(y,t)G s (x - y)dy. 



_* 2 

In this paper, we use a Gaussian filter as in [3], £,5(3;) = . 
On convolving ([I]) with Gs, the large eddy solution u is to satisfy 

— /"* 1 _ 

«t = u xx + u - u 3 + ift u ( x > s )ds, 



1 + It - s|0 



or 



_„ /" f 

«t = u xx + u-u ,i + I -g u(x,s)ds + R(x,t), (2) 

Jo 1 + |t — s| p 

where the remainder term, i.e., the subgrid scale (SGS) term R(x,t) is defined as 



R(x,t) := (u) 6 - (u 3 ). (3) 

We can write u = u + u' with u the large eddy term and u' the fluctuating term. Note 
that u = u — v! . So the SGS term R(x,t) involves nonlinear interactions of fluctuations 
v! and the large eddy flows. Thus R(x,t) may be regarded as a function of u and u': 
R := R(u,u'). 

The leads to a possibility of approximating R(x, t) by a suitable stochastic process 
defined on a probability space (0,JF, P), with u 6 f2, the sample space, cr— field ^ and 
probability measure P. This means that we treat R data as random data as in [22], which 
take different realizations, e.g., due to fluctuating observations or due to numerical simu- 
lation with initial and boundary conditions with small fluctuations. In fluid or geophysical 
fluid simulations, the SGS term may be highly fluctuating and time-correlated [22], and 
this term may be inferred from observational data \27\ 128]. or from fine mesh simulations. 



3 



In fact, in our case study here, the subgrid scale term R(x, t) is clearly time-correlated; 
see Fig. 1. The (averaged) time correlation function here, over a computational time 
interval [0, T], is defined as : 



Corr(x, s) 



cov(R(x, t), R(x, t + s)) 



-dt, 



T J STD{R(x,t)) ■ STD(R(x,t + s)) 
where cov denotes covariance, and STD(R) = \/E(i? — Ei?) 2 is the standard deviation. 



1.02 




Figure 1: Corr(0,s) — Averaged time correlation of the subgrid scale term R(x,t), at x = 0, 
for u t = u xx + u — u 3 + J ppr^jra u(x, s)ds, u(x,Q) = 0.53a; — 0.47sin(1.57ra;), u(—l,t) = 
-1, u(M) = l;/3 = 2 

This further suggests for parameterizing the subgrid scale term R(x, t) as a time- 
correlated or colored noisy term. 

Before we devise how we parameterize the subgrid scale term R(x, t) as a colored noise 
or time-correlated term, we first discuss fractional Brownian motion and colored noise in 
the next section. 



3 Stochastic parameterization via a colored noise 

We first discuss a model of colored noise in terms of fractional Brownian motion. The 
fractional Brownian motion is a generalization of the more well-known process of Brown- 
ian motion. It is a centered Gaussian process with stationary increments. However, the 
increments of the fractional Brownian motion are not independent, except in the standard 
Brownian case (H = ^). The dependence structure of the increments is modeled by a so 
called Hurst parameter H S (0, 1). For more details, see [2T | [23 1 [6l [TBI [35] . 
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Definition of fractional Brownian motion: For H G (0, 1), a Gaussian process B H (t), 
or fBM(t), is a fractional Brownian motion if it starts at zero B H (0) = 0, a.s., has mean 
zero E[B H (t)] = 0, and has covariance E[B H (t)B H {s)} = ±(\t\ 2H + \s\ 2H - \t - s\ 2H ) for 
all t and s. The standard Brownian motion is a fractional Brownian motion with Hurst 
parameter H = ^. 

Some properties of fractional Brownian motion: A fractional Brownian motion B H (t) 
has the following properties: 

(i) It has stationary increments; 

(ii) When H = 1/2, it has independent increments; 

(hi) When H ^ 1/2, it is neither Markovian, nor a semimartingale. 

The exact simulation of B H (t\), .., B H (t n ) is in general computationally very expensive. 
The Cholesky decomposition method, which is to our knowledge a known exact method 
for the non-equidistant simulation of fractional Brownian motion, requires 0(n 3 ) opera- 
tions. Moreover the covariance matrix, which has to be decomposed, is ill-conditioned. 
If the discretization is equidistant, i.e. , t{ = i/n,i = l,...,n, the computational cost 
can be lowered considerably. For example, the Davis-Harte algorithm for the equidistant 
simulation of fractional Brownian motion has computational cost 0(nlog(n)); see, e.g., 
Craigmile [I]. 

Here we use the Weierstrass-Mandelbrot function to approximate the fractional Brow- 
nian motion. The basic idea is to simulate fractional Brownian motion by randomizing 
a representation due to Weierstrass. Given the Hurst parameter H with < H < 1, we 
define the function w(t) to approximate the fractional Brownian motion: 



where r = 0.9 is a constant, Cj's are normally distributed random variables with mean 
and standard deviation 1, and the dj's are uniformly distributed random variables in the 
interval < dj < 2ir. The underlying theoretical foundation for this approximation can 
be found in |29^ I20j. Fig. 2 shows a sample path of the fractional Brownian motion, when 



The increments of fractional Brownian motion are correlated in time. This motivates us 
to parameterize the subgrid scale term R(x,t), which is time-correlated, by using colored 
noise B^ . 

We thus parameterize the subgrid scale term R(x, t) term ([3]) above as a mean com- 
ponent plus a colored noise. To be more specific, we model R(x,t) as follows: 



oc 




Hurst parameter H = |. 



R(x,t) = f{u)+cr{x) 



dB? 



(4) 



dt 



where 



dt 



is a colored noise, and 



f[u) =ER(x,t) 



(5) 
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Figure 2: A sample path of fractional Brownian motion B 



(t), with H = 0.75 



is the mean component of the subgrid scale term R(x,t). Moreover, the noise intensity 
cr(x) is a non-negative deterministic function to be determined from fluctuating SGS data 
R. The subgrid scale term R(x, t) may be inferred from observational data [27} I28|. or from 
fine mesh simulations as we do here. We represent the mean component f(u) in terms of 
the large eddy solution u. The specific form for / depends on the nature of the mean of 
R. Here we take f(u) = ao + a\u + 02U 2 + a^u 3 , where coefficients a,'s are determined via 
data fitting by minimizing J D [ao + a\u + (12U 2 + a^u 3 — KR(x, t)] 2 dxdt. Moreover, we 
take Bf as a scalar fractional Brownian motion. 

Note that a is to be calculated or estimated from the fluctuating SGS data R, either 
from observation or (in this paper) from fine mesh simulations; see detailed discussions in 
[221 [5]. So this is an inverse problem. As in usual inverse problems [33], the stochastic 
parameterizations for the SGS term R is not unique. This offers an opportunity for 
trying various stochastic parameterization schemes, much as one uses various smoother 
functions (e.g., polynomials or Fourier series) to approximate less regular functions or data 
in deterministic approximation theory. 

To estimate the unknown parameter (function) cr(x), we start with @-([5]) to get the 
following relation: 



Taking time integral over a computational interval [0, T] on both sides, we obtain 



R(x,t) -ER(x,t) 




dBf 



(6) 



dt 



/ [R(x, t) - ER(x, t)]dt = / a(x)dBf I = a{x)B% 



Jo Jo 



Therefore, taking mean-square on both sides, 



E( [ [R(x,t) -ER(x,t)]dt) 2 = a 2 {x)T 2H . 
Jo 

Thus an estimator for cr(x) is 



a(x) = -LjE(£[R(x,t) -ER(x,t)]dt)i , (7) 

which can be computed numerically. 

By the stochastic parameterization Q on the SGS term R, with / determined from 
© and a from ([7]), the LES model ([2]) becomes a stochastic partial differential equation 
(SPDE) for the large eddy solution U ~ u: 

rt i dB H 

U t = U xx + U -U' 3 + / — — rg U(x, s)ds + f(U) + C7(X)—^-, (8) 

Jo 1 + |i — s\^ at 

with boundary conditions U{— 1, t) = a, U(l,t) = b and filtered initial condition 

U(x,0)=u (x). (9) 



4 Numerical Experiments 

We use a spectral method to solve nonlinear system ([T]) and (JH]) numerically. For more 
details, please see [M]. We take the following initial and boundary conditions: 

u(x, 0) = uq = 0.53a; — 0.47sm(1.57ra;), u(—l,t) = —l, u(l,t) = 1 

Fine mesh simulations of the original system with memory (TjQ) are conducted to gener- 
ate benchmark solutions or solution realizations, with initial conditions slightly perturbed; 
see Fig. 3. These fine mesh solutions u are used to generate the SGS term R defined in 
(J3J) at each time and space step. The filter size used in calculating R is taken as 5 = 0.01. 
The mean / is calculated from ([5]) via cubic polynomial data fitting (as discussed in the 
last section), and parameter function cr(x) is calculated as in ([7]). The stochastic LES 
model (|8|) is solved by the same numerical code but on a coarser mesh. Note that a four 
times coarser mesh simulation with no stochastic parameterization for the original system 
(H|) does not generate satisfactory results; see Fig. 4. The stochastic LES model (jSJ) is 
then solved in the mesh four times coarser than the fine mesh used to solve the original 
equation ([T]). The stochastic parameterization leads to better resolution of the solution 
as shown in Fig. 5. The root-mean-square error, error(x,t) := ^~E\u(x,t) — U(x,t)\ 2 , 
between the fine mesh solution u (Fig. 3) and this stochastic LES solution U (Fig. 5) is 
plotted in Fig. 6. 
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Figure 3: Solution to the original system on a fine mesh, u t = u xx + u — u 3 + J * 1+ ^_ s ^ u(x, s)ds, 
(3 = 2, mesh size Ax = 0.001. 




Figure 4: Solution to the original system with NO stochastic parametrization on the mesh four 
times coarser than the mesh used in Fig. 3, u t = u xx + u — u 3 + f Q 1+ i^_ s w u(x, s)ds, (3 = 2. 
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Figure 5: Solution to LES model with stochastic parametrization on the mesh four times coarser 
than the mesh used in Fig. 3, U t = U XX + U - U 3 + f* „ U(x, s)ds + f(U) + a{x)B^ , [3 = 2, 




Figure 6: Mean-square error of the LES model on the mesh four times coarser than the mesh used 
in Fig. 3, U t = U xx + U-U a + f* * ,„ U{x, s)ds + f(U) + a{x)B» , /3 = 2, H = §. 
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5 Discussions 



We have discussed the issue of modeling the impact of unresolved scales on the resolved 
ones, in the context of large eddy simulation of a nonlinear partial differential equation 
with memory. The resulting model is a stochastic partial differential equation, which 
describes large scale evolution with some effects of unresolved scales taken into account. 

It would be very interesting to investigate whether this stochastic approach works 
for simulation of other complex systems, such as climate systems, fluid flows, biological 
systems and materials. 
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